Sei’s drum!
Michael Großbach
Draft 2020-08-21
library(brms)
library(dplyr)
library(ggplot2)
library(ggthemes)
library(kfigr)
library(knitr)
library(patchwork)The data set was read in as is, the columns subj, arm, cond, series and no were converted to factors, the column filename was deleted. A new column was added to the dataset, patient, based on the subj prefix (‘P’ or ‘S’), with a ‘0’ coding for healthy controls, and a ‘1’ for patients.
read.csv("features.csv") %>%
mutate(filename = NULL,
patient = factor(if_else(grepl("^P",
subj),
1L,
0L)),
series = factor(series),
subj = factor(subj),
arm = factor(arm),
cond = factor(cond),
no = factor(no)) ->
drum_beats1 Exploration
1.1 Instruction Condition, and Arm Used
In a first attempt to get an overview, the descriptors’ density estimates were plotted, broken down by arm and cond.
dep_vars <- names(drum_beats)[which(!(names(drum_beats) %in% c("subj",
"arm",
"cond",
"no",
"series",
"patient")))]
n_vars <- length(dep_vars)
plot_lst <- vector("list", length = n_vars)
plt_cnt <- 1
for (dv in dep_vars) {
if (plt_cnt == n_vars) {
plt_pos <- c(2, 0.5)
} else {
plt_pos <- "none"
}
if (plt_cnt %% 4 == 1) {
ylab <- "density"
} else {
ylab <- ""
}
plot_lst[[plt_cnt]] <- ggplot(drum_beats,
aes_string(x = dv,
color = "cond")) +
geom_density(na.rm = TRUE) +
scale_color_colorblind() +
ylab(ylab) +
geom_rug(alpha = 0.25) +
theme(legend.position = plt_pos) +
facet_wrap(~ arm)
plt_cnt <- plt_cnt + 1
}
print(wrap_plots(plot_lst,
ncol = 4) +
plot_annotation(
tag_levels = "A"))
Fig. 1. Descriptor densities along with raw data points (rug ticks on x axes), broken down by arm (dominant [D] vs non-dominant [ND]) and instruction condition (controlled = C, normal = N).
Several peculiarities are immediately apparent from the plots in Fig. 1:
- There are no substantial differences in the density distributions between normal (N) and controlled (C) strokes
attDur,LAT, andattFlathave very few unique data values (see rug ticks at bottom of plots)- Looking at the dominant arm alone
(and excluding the vars, the densities have a trend toward a larger variance in the controlled condition (in other words: a less pointed distribution) with additional or larger humps in the tails, respectivelyattDur,LAT, andattFlatfor the reason just mentioned)
The lack of any clear-cut differences in the raw data between conditions suggests that it will be hard to find any differences by modeling.
The limited number of unique values in the attack-related measures (e.g. attDur only has 35 unique values in a total of 1102 observations across all subjects, conditions, and sides; that’s only 3 percent!) suggests that rounding errors propagated through the calculations. This is probably due to very similar, i.e. highly automated, and short attack times combined with the given sampling frequency, resulting in few data points, which in course of the calculations result in the observed phenomenon. But I’m just guessing here, as I do not know anything about these descriptors and how they are calculated or interpreted. Judging by the range of values of, i.e., attDur (2.59, 9.8ms, across both sides and conditions) and a sampling frequency of 48 kHz, this leaves us at 48 * (9.8 - 2.59) = 346 points to choose beginning and end of the attack. Given the supposed highly automatized motor program used to initiate the stroke, along with the laws of physics at play here (no pun intended), it is not very surprising to see very few unique values. My limited knowledge—or rather, my ignorance of anything sound-related—aside, from a statistical standpoint these measures do not seem suited to describe any differences between the experimental conditions investigated here.
Given the very few data points in the attack phase, any descriptor derived from such a short period of time (.e.g.. attSPL) cannot be judged as being stable in the sense of being reproducible. Hence I suggest to drop all attack-derived descriptors.
# drum_beats %>%
# mutate(attDur = NULL,
# LAT = NULL,
# attSPL = NULL,
# attSC = NULL,
# attFlat = NULL,
# attSpecFlat = NULL) ->
# drum_beatsThe increase in variance for the conditions N < C does not come as a surprise as I assume normal also means highly trained and thus automatized, whereas controlled involves less automatization and more ‘individualness’ both within and between subjects.
1.1.1 Conclusions
Looking at the above investigated descriptors and Danielsen et al. (2015), it seems reasonable to limit the modeling attempts to totDur, totSPL, totSC, and TC.
But Sofia wrote on 2020-07-09: “Francesco and I have discussed a bit related to descriptors and we want to concentrate on the “transient” period (although the name probably will change). Spectral Centroid (transSC) should be one, and I suggest transFlat for the other. Francesco, does that sound reasonable?"
On 2020-07-22 both Sofia and Francesco agreed upon transSC, transFlat, and transCrest as the probably most important response variables to look at.
1.1.2 Francesco’s comments
Update:
- The lack of data points for
attFlatis fixed by solving a bug in the feature extraction: now the downsampling ratio for the envelope extraction is reduced, and makes the MIRToolbox algorithm less sensitive to frames with an rms value of 0 (which returns an incorrect value of 0, since we have a geometric mean at the numerator). Now we have a larger spread which allows high values (1 = peaky envelope, 0 = smooth/flat envelope). - The crest factor is now calculated for the separate phases.
Regarding the attack phase: I agree with the remarks regarding duration. Given the fact that we are not comparing different instruments, I wouldn’t have expected a large variation. This is not surprising if we consider that our system is changing only slightly (same rototom, same drumstick, same action, a bit different tuning across subjects): in fact, the perceived timbral differences are so small that we are in trouble guessing on the descriptors.
The sampling frequency is even lower (\(f_s = 44100\) Hz). Even if we had a higher sample rate which could reveal some discrepancies in the attack durations, we would have to prove that they are perceptually relevant.
Therefore, I am happy to discard attDur and LAT (which is obviously a log-transformed duration, only there for the sake of consistency with the literature).
I am a bit more in doubt when it comes to discarding all the attack descriptors. Even if what Michael says is true from a statistical point of view, we should still be able to discriminate timbre on short time windows due to the high temporal resolution of our hearing. Attack phase descriptors (with the same definition of attack that we are using, which is most likely not coincident with perceptual attack) are employed in Câmara et al. (2020), and the Oslo group has a paper under construction which analyzes drum sounds in a similar manner (see OsloPlots.png in the OneDrive folder).
I am worried that merely taking the overall descriptors into account would introduce a lot of unnecessary and perceptually catching information — mostly the tonal part of the signal, i.e. the drum ringing in the last part of the decay. That’s why Sofia and I are suggesting to look at what we could call “main energy” or “early decay” phase (i.e. from max peak to temporal centroid).
Would it be feasible to set up 4 different models (i.e. one for each phase), at least in the univariate version?
As for the descriptors to include: although TC is employed in Danielsen et al. (2015), the PDFs are even more similar. I would go for Dur (except attack?), SPL, SC, and one between Flat, SpecFlat or Crest.
My (informal and biased) listening tests tell me that, at least for some subjects, I hear a pattern going towards a harsher (controlled) vs smoother (normal) timbre exactly at the hit point, plus slightly less (controlled) or more pitch/amplitude fluctuation. Hopefully this could be catched by spectral centroid, specrtral/temporal flatness, or crest factor. This should be independent of SPL unless the subject misinterpreted the instructions, therefore SPL acts as a sort of control variable in our model.
1.1.3 Model Considerations
There are several points that need to be considered before making a decision regarding the type of modeling to be done in this study, (1) the sample size places restrictions on the external validity; (2) data from small samples can be better modeled when regularization is in place to ‘tame’ the estimates; (3) the hierarchical structure of the data (subjects played several trials with either their dominant or non-dominant arm under two conditions) suggests a multilevel analysis of the data which would, in addition to the Bayesian regularization via priors, also results in shrinkage of the estimates; (4) given the small distribution differences between the two experimental conditions, along with the
- Given the extremely small sample size, any attempt to fit the data using traditional statistics, a.k.a. null-hypothesis significance testing, would make us vulnerable to all kinds of criticism. I therefore suggest to use Bayesian statistics instead, as it allows us to estimate probability distributions of parameters rather than confidence intervals around point estimates, and thus embraces uncertainty in estimates.
- Additionally, Bayesian regression uses prior probability distributions to arrive at sensible estimates. These priors regularize estimates, or draw them toward zero, a desirable effect which has long been recognized even in traditional statistics (e.g. ridge regression and lasso; Tibshirani (1996)).
- ANOVA (or regression with grouping variables) requires each subject’s trials to be averaged, say, within conditions or groups (or both), to be able to assess treatment effects; this results in the individual variation within a subject to get lost while it might have been very informative to include it in the analysis. Multilevel (or hierarchical) modeling allows to include the entire data structure (Fig. 2) so that no information gets lost through averaging, but all variation (both individual and treatment-driven) propagates through the analysis and ends up in the final results, allowing for more realistic credibility margins around estimates.
- Variance in traditional statistics is considered to be fixed (think homoscedasticity in ANOVA, or the \(\epsilon\) in regression formulas such as \(y_i = \beta_0 + \beta_1x_i + \epsilon\)), whereas it is an estimated quantity in Bayesian statistics and therefore can vary between conditions (or groups or subjects etc.) and thus capture varying variance in different conditions, or groups etc.
# chunk intentionally empty Group: Patient Healthy Control
/ \ / \
Instruction: Normal Controlled Normal Controlled
| | \ / | \
Player: 1 2 ... n 1 2 ... n
/ \
Side dominant non-dominant
/ | \ / | \
Trial 1 .. p 1 .. pFig. 2. The multilevel structure of a data set should be reflected in the analysis.
1.1.4 Conclusions
Will start with a simple univariate model, add predictors and interactions, then a bivariate model, and finally a quadruple-variate model and see where this leads us.
1.2 Individual Responses
The humps in the density distributions in section Instruction Condition, and Arm Used made me curious where they might originate from. So in the following graph the density estimates of one of the descriptors, transSC, are plotted broken down by subj. and then also by cond, separately for patients and healthy subjects.
ggplot(drum_beats, aes(transSC,
color = subj,
linetype = cond)) +
geom_density(alpha = 0.1) 
Fig. 3. transSC density distribution of participants, broken down by instruction condition.
It is obvious that some participants do not differ substantially between the normal and the controlled condition, whereas others do, and even markedly so. Additionally, there seems to be quite a spread of the centers of distributions across a wide range of values, suggesting very individual drum sounds.
The wide distribution of centers of mass between individuals made me want to further break down the plot.
ggplot(drum_beats, aes(transSC,
color = subj,
linetype = cond)) +
geom_density(alpha = 0.1) +
facet_wrap(~ patient)
Fig. 4. transSC density distribution of participants, broken down by instruction condition and group (patients and healthy subjects).
The above plot is interesting in that it seems to show that the four healthy subjects were more uniform than the patients in their transSC distributions and also, with the exception of S2, had very similar transSC distributions for the normal and the controlled conditions. In the patients, two had very similar distributions in both conditions (P1 and P5), whereas the two others showed differing results for the two conditions. So with regard to Sofia’s hypothesis (ch. Modeling) I’d argue, at least for transSC alone, having drummers play normal and controlled strokes would not allow subjects to tell the difference in a listening test. But maybe it would qualify as a screening test for movement disorders in drummers. Just a random thought.
While breaking it down it occurred to me that looking at individual variation (that is, between series) might also be enlightening.
ggplot(drum_beats, aes(transSC,
color = cond,
group = series)) +
geom_density(alpha = 0.1) +
facet_wrap(~ subj, nrow = 2)
Fig. 5. Density plots of individual series.
And it was! The plot in Fig. 5 made visible that P1 and P5 had consistently very low transSC values. P3 was consistent within the normal instruction condition, but had, on average, higher values, and with a lot more variation, in the controlled condition. P4 had more variation in both conditions, and higher transSC values in both; in other words, P4 was consistantly bad. (But then again: what do I know what bad is wrt transSC!)
The healthy subjects showed comparable variation and centers of mass within conditions, and all but one (S2) also across conditions.
Comparing the two rows of panels in Fig. 5 reveals that healthy subjects have more variation than the patients.
Although tempting, we probably shouldn’t get carried away and generalize to the populations of healthy drummers and ones with movement disorders, respectively.
2 Modeling
Sofia, on 2020-07-07: “The main hypothesis is that playing instruction (N/C) will affect the stroke in a way that is perceivable”.
2.1 Univariate Models
Looking at the variables agreed upon I have decided to use as response variables in the regression models—implying that this is by no means set in stone—, it seems like a skewed normal link function would be appropriate to model them.
Let’s start with transSCtotDur
(m0_form <-bf(transSC ~ 1 + (1 | subj)))transSC ~ 1 + (1 | subj) which claims that transSC is explained by (‘~’) an intercept, denoted by ‘1’, and an additional term ‘(1 | subj)’. The ‘1’ in parentheses again stands for the intercept, but the pipe ‘|’ assigns an intercept to each level of the factor ‘subj’. In this particular case this means that the model will estimate an individual intercept for each unique drummer listed in the data set column subj. These individual, or varying, intercepts are then used in informing the estimation of the population intercept.
Note: set MODEL to TRUE at the top of the script if you didn’t compile/build your model yet.
if (MODEL) {
m0 <- brm(m0_form,
family = skew_normal(),
inits = "0",
data = drum_beats)
m0 <- add_criterion(m0,
"loo",
reloo = TRUE)
save(m0,
file = "m0.rda")
} else {
load("m0.rda")
}This null model is also termed unconditional model because it has no grouping structure apart from individuals–the ‘subj’ bit in the model equation above. There is some variation in every natural data set. To make sure, it’s not just variation caused by different participants, we can calculate the intra-class correlation coefficient (ICC).
m0_icc <- ICC(m0, "subj")The Null model’s ICC amounts to 0.75, which suggests that approx. 75 percent of the variation in the data set can be attributed to (or explained by) the grouping structure. This is highly unfortunate, as it does not leave a lot of variation to be explained by independent factors like instruction, or arm. The most likely reason for this high ICC value is the small sample size combined with high inter-individual variation. Small sample sizes combined with large trial numbers are less of a problem when subjects respond, on average, close to the population mean, even with large spread due to fluctuating alertness, increasing fatigue etc., .e.g. in reaction time paradigms. But here, with large intra- and inter-individual variation, this might become a problem.
Tab. 1. Model summary.
(m0_summary <- summary(m0,
priors = TRUE)) Family: skew_normal
Links: mu = identity; sigma = identity; alpha = identity
Formula: transSC ~ 1 + (1 | subj)
Data: drum_beats (Number of observations: 1102)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
Priors:
alpha ~ normal(0, 4)
Intercept ~ student_t(3, 735, 136.5)
sd ~ student_t(3, 0, 136.5)
sigma ~ student_t(3, 0, 136.5)
Group-Level Effects:
~subj (Number of levels: 8)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 113.47 32.93 68.06 193.47 1.01 674 1254
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 746.78 37.92 670.61 821.58 1.01 717 1226
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 64.94 1.57 61.84 68.01 1.00 2080 1822
alpha 3.14 0.35 2.49 3.84 1.00 1930 1791
Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).The intercept in Tab. 1 amounts to roughly 747. In this model specification, the intercept is identical to the data set average. The table also shows that the inter-individual standard deviation (sd(Intercept)) is large compared to the unexplained variation (sigma). This led to the large ICC value above. By adding more and more independent factors to the model specification, we will later try to decrease \(\sigma\), i.e. ‘explain away’ as much remaining variation as possible.
pp_check(m0, nsamples = 100)
Fig. 6. Posterior predictive check. The thick blue line shows the distribution of the empirical data. The thin blue lines are one-hundred realizations of data generated from parameters estimated by the model.
The posterior predictive plot in Fig. 6 gives an impression on how the null model would generate data, given the parameters it estimated from the empirical data. As the thick line deviates from the modeled thin lines, especially on the left and the right side of the peak of the distribution, it is apparent that the model can be improved.
# conditional_effects(m0) 2.1.1 Group Model
In this model, the intercept is complemented with second ‘main’ or population effect, the grouping variable patient:
(m1_form <- bf(transSC ~ 1 + patient +
(1 | subj)))transSC ~ 1 + patient + (1 | subj) I could have also specified the model without the explicit ‘1 +’, as the intercept is implicitly included in the model unless I explicitly exclude it. From now on I will always save some extra typing by refraining from explicitly indluding the intercept in models.
So now the model not only contains the individuals as grouping structure to ‘explain away’ variation, but also whether they belong to the patients or the healthy subjects.
The prior distributions in Bayesian models reflect the knowledge about the estimated parameters. The package brms automatically places weakly informative priors on parameters as soft contraints. But with more complex models involving varying parameter estimates, the statistical back end which does the heavy lifting, needs stronger priors particularly on these parameters, otherwise models don’t converge. The parameters most vulnerable to outliers in the data are the varying effects parameters, or random effects in traditional statistics. Therefore we place a stronger prior probability distribution over the estimate of the SD of individual intercepts:
(m1_prior <- set_prior("normal(0, 10)", class = "sd"))sd ~ normal(0, 10)and leave the rest of the priors as suggested by the package brms (see Tab. 2 for their priors).
if (MODEL) {
m1 <- brm(m1_form,
prior = m1_prior,
inits = "0",
family = skew_normal(),
data = drum_beats)
m1 <- add_criterion(m1,
"loo",
reloo = TRUE)
save(m1,
file = "m1.rda")
} else {
load("m1.rda")
}Tab. 2. Model summary.
(m1_summary <- summary(m1,
priors = TRUE)) Family: skew_normal
Links: mu = identity; sigma = identity; alpha = identity
Formula: transSC ~ 1 + patient + (1 | subj)
Data: drum_beats (Number of observations: 1102)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
Priors:
alpha ~ normal(0, 4)
Intercept ~ student_t(3, 735, 136.5)
sd ~ normal(0, 10)
sigma ~ student_t(3, 0, 136.5)
Group-Level Effects:
~subj (Number of levels: 8)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 49.08 4.87 40.10 59.23 1.00 2584 2541
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 761.13 25.06 710.46 810.22 1.00 1130 1658
patient1 -23.30 35.05 -92.84 47.46 1.00 1255 1588
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 65.10 1.54 62.07 68.15 1.00 3728 2679
alpha 3.26 0.34 2.63 3.96 1.00 3904 2723
Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).The intercept in Tab. 2 amounts to roughly 761. In this model specification, the intercept is identical to the controlled strokes, while the patient1 value is the mean of the posterior distribution difference between the healthy subjects and the patients (-23). The table also shows that the interindividual standard deviation (sd(Intercept) \(\approx\) 49 is still large compared to the unexplained variation (sigma \(\approx\) 65).
pp_check(m1, nsamples = 100)
Fig. 7. Posterior predictive check. The thick blue line shows the distribution of the empirical data. The thin blue lines are one-hundred realizations of data generated from parameters estimated by the model.
The posterior predictive plot in Fig. 7 gives an impression on how this model would generate data, given the parameters it estimated from the empirical data. There is not much change from Fig. 6, which is not very surprising given the small estimated difference between the groups.
conditional_effects(m1)
2.1.2 Group + Condition Model
In this model, the intercept is complemented with second ‘main’ or population effect, the instruction condition cond:
(m2_form <-bf(transSC ~ patient + cond +
(1 | subj)))transSC ~ patient + cond + (1 | subj) So now the model also contains the individuals as grouping structure to ‘explain away’ variation, but also the manipulations.
The prior distributions in Bayesian models reflect the knowledge about the estimated parameters. The package brms automatically places weakly informative priors on parameters as soft contraints. But with more complex models involving varying parameter estimates, the statistical back end which does the heavy lifting, needs stronger priors particularly on these parameters, otherwise models don’t converge. The parameters most vulnerable to outliers in the data are the varying effects parameters, or random effects in traditional statistics. Therefore we place a stronger prior probability distribution over the estimate of the SD of individual intercepts:
(m2_prior <- c(set_prior("normal(0, 3)", class = "sd"),
set_prior("normal(0, 3)", class = "sigma"),
set_prior("normal(0, 2)", class = "alpha"))) prior class coef group resp dpar nlpar bound source
normal(0, 3) sd user
normal(0, 3) sigma user
normal(0, 2) alpha userand leave the rest of the priors as suggested by the package brms (see Tab. 3 for their priors).
if (MODEL) {
m2 <- brm(m2_form,
prior = m2_prior,
family = skew_normal(),
inits = "0",
data = drum_beats)
m2 <- add_criterion(m2,
"loo",
reloo = TRUE)
save(m2,
file = "m2.rda")
} else {
load("m2.rda")
}Tab. 3. Model summary.
(m2_summary <- summary(m2,
priors = TRUE)) Family: skew_normal
Links: mu = identity; sigma = identity; alpha = identity
Formula: transSC ~ patient + cond + (1 | subj)
Data: drum_beats (Number of observations: 1102)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
Priors:
alpha ~ normal(0, 2)
Intercept ~ student_t(3, 735, 136.5)
sd ~ normal(0, 3)
sigma ~ normal(0, 3)
Group-Level Effects:
~subj (Number of levels: 8)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 29.28 1.53 26.40 32.46 1.00 2941 2444
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 780.34 14.97 751.58 811.34 1.00 1500 1967
patient1 -27.76 21.28 -70.57 13.29 1.00 1621 1999
condN -41.31 3.39 -47.91 -34.88 1.00 3868 3201
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 52.36 0.94 50.61 54.18 1.00 4011 2896
alpha 1.79 0.22 1.37 2.25 1.00 3857 3095
Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).The intercept in Tab. 3 amounts to roughly 780. In this model specification, the intercept is identical to the controlled strokes, while the condN value is the mean of the posterior distribution difference between controlled and normal strokes (-41). The table also shows that the interindividual standard deviation (sd(Intercept) \(\approx\) 29 is still large compared to the unexplained variation (sigma \(\approx\) 52).
pp_check(m2, nsamples = 100)
Fig. 8. Posterior predictive check. The thick blue line shows the distribution of the empirical data. The thin blue lines are one-hundred realizations of data generated from parameters estimated by the model.
The posterior predictive plot in Fig. 8 gives an impression on how this model would generate data, given the parameters it estimated from the empirical data. There is not much change from Fig. 6, which is not very surprising given the small estimated difference between the conditions.
conditional_effects(m2)
2.1.3 Group + ConditionVI Model
The last model included the instruction condition, more realistically reflecting the true structure of the data set. But it did not acknowledge that each subject executed several strokes in each of these conditions. This model includes a term with a varying intercept (VI) for condition to reflect just that, which will also assist in more realistically estimate the population effect of cond:
(m3_form <- bf(transSC ~ patient + cond +
(1 | subj) +
(1 | cond)))transSC ~ patient + cond + (1 | subj) + (1 | cond) Including more varying parameters in the model requires the prior on them to be even stronger:
(m3_prior <- m2_prior) prior class coef group resp dpar nlpar bound source
normal(0, 3) sd user
normal(0, 3) sigma user
normal(0, 2) alpha userif (MODEL) {
m3 <- brm(m3_form,
prior = m3_prior,
family = skew_normal(),
inits = "0",
data = drum_beats)
m3 <- add_criterion(m3,
"loo",
reloo = TRUE)
save(m3,
file = "m3.rda")
} else {
load("m3.rda")
}Tab. 4. Model summary.
(m3_summary <- summary(m3,
priors = TRUE)) Family: skew_normal
Links: mu = identity; sigma = identity; alpha = identity
Formula: transSC ~ patient + cond + (1 | subj) + (1 | cond)
Data: drum_beats (Number of observations: 1102)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
Priors:
alpha ~ normal(0, 2)
Intercept ~ student_t(3, 735, 136.5)
sd ~ normal(0, 3)
sigma ~ normal(0, 3)
Group-Level Effects:
~cond (Number of levels: 2)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 2.39 1.77 0.10 6.59 1.00 2228 1784
~subj (Number of levels: 8)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 29.28 1.51 26.43 32.30 1.00 2915 2550
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 779.70 15.39 749.79 810.03 1.00 1429 2093
patient1 -27.10 21.30 -68.53 14.14 1.00 1817 2100
condN -41.22 5.58 -52.39 -30.63 1.00 2578 2255
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 52.37 0.94 50.59 54.20 1.00 4592 2856
alpha 1.80 0.23 1.36 2.27 1.00 3963 2284
Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).The intercept in Tab. 4 amounts to roughly 780. In this model specification, the intercept is identical to the controlled strokes, while the condN value is the mean of the posterior distribution difference between controlled and normal strokes (-41). The table also shows that the inter-individual standard deviation (sd(Intercept) \(\approx\) 29 is not large anymore compared to the unexplained variation (sigma \(\approx\) 52).
pp_check(m3, nsamples = 100)
Fig. 9. Posterior predictive check. The thick blue line shows the distribution of the empirical data. The thin blue lines are one-hundred realizations of data generated from parameters estimated by the model.
The posterior predictive plot in Fig. 9 gives an impression on how this model would generate data, given the parameters it estimated from the empirical data. There is not much change from Fig. 6, which is not very surprising given the small estimated difference between the conditions.
conditional_effects(m3)
2.1.4 Group x ConditionVI
The last model included the instruction condition, more realistically reflecting the true structure of the data set. But it did not acknowledge that each subject executed several strokes in each of these conditions. This model includes a term with a varying intercept (VI) for condition to reflect just that, which will also assist in more realistically estimate the population effect of cond:
(m4_form <- bf(transSC ~ patient * cond +
(1 | subj) +
(1 | cond)))transSC ~ patient * cond + (1 | subj) + (1 | cond) Including more varying parameters in the model requires the prior on them to be even stronger:
(m4_prior <- m3_prior) prior class coef group resp dpar nlpar bound source
normal(0, 3) sd user
normal(0, 3) sigma user
normal(0, 2) alpha userif (MODEL) {
m4 <- brm(m4_form,
prior = m4_prior,
family = skew_normal(),
inits = "0",
data = drum_beats)
m4 <- add_criterion(m4,
"loo",
reloo = TRUE)
save(m4,
file = "m4.rda")
} else {
load("m4.rda")
}Tab. 5. Model summary.
(m4_summary <- summary(m4,
priors = TRUE)) Family: skew_normal
Links: mu = identity; sigma = identity; alpha = identity
Formula: transSC ~ patient * cond + (1 | subj) + (1 | cond)
Data: drum_beats (Number of observations: 1102)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
Priors:
alpha ~ normal(0, 2)
Intercept ~ student_t(3, 735, 136.5)
sd ~ normal(0, 3)
sigma ~ normal(0, 3)
Group-Level Effects:
~cond (Number of levels: 2)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 2.39 1.83 0.07 6.88 1.00 2373 1368
~subj (Number of levels: 8)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 29.25 1.52 26.35 32.38 1.00 2682 2525
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 780.52 15.72 750.05 812.34 1.00 1723 1958
patient1 -27.98 21.77 -69.25 14.30 1.00 1703 2297
condN -42.29 6.13 -54.39 -29.83 1.00 2628 1939
patient1:condN 2.24 6.37 -10.51 14.61 1.00 3979 2902
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 52.37 0.92 50.62 54.28 1.00 4866 3208
alpha 1.81 0.23 1.35 2.27 1.00 4122 2687
Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).The intercept in Tab. 5 amounts to roughly 781. In this model specification, the intercept is identical to the controlled strokes, while the condN value is the mean of the posterior distribution difference between controlled and normal strokes (-42). The table also shows that the inter-individual standard deviation (sd(Intercept) \(\approx\) 29 is not large anymore compared to the unexplained variation (sigma \(\approx\) 52).
pp_check(m4, nsamples = 100)
Fig. 10. Posterior predictive check. The thick blue line shows the distribution of the empirical data. The thin blue lines are one-hundred realizations of data generated from parameters estimated by the model.
The posterior predictive plot in Fig. 10 gives an impression on how this model would generate data, given the parameters it estimated from the empirical data. There is not much change from Fig. 6, which is not very surprising given the small estimated difference between the conditions.
conditional_effects(m4)
2.1.5 Group x ConditionVI + Arm Model
(m5_form <- bf(transSC ~ patient * cond + arm +
(1 | subj) +
(1 | cond)))transSC ~ patient * cond + arm + (1 | subj) + (1 | cond) (m5_prior <- m4_prior) prior class coef group resp dpar nlpar bound source
normal(0, 3) sd user
normal(0, 3) sigma user
normal(0, 2) alpha userif (MODEL) {
m5 <- brm(m5_form,
prior = m5_prior,
family = skew_normal(),
inits = "0",
data = drum_beats)
m5 <- add_criterion(m5,
"loo",
reloo = TRUE)
save(m5,
file = "m5.rda")
} else {
load("m5.rda")
}Tab. 6. Model summary.
(m5_summary <- summary(m5,
priors = TRUE)) Family: skew_normal
Links: mu = identity; sigma = identity; alpha = identity
Formula: transSC ~ patient * cond + arm + (1 | subj) + (1 | cond)
Data: drum_beats (Number of observations: 1102)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
Priors:
alpha ~ normal(0, 2)
Intercept ~ student_t(3, 735, 136.5)
sd ~ normal(0, 3)
sigma ~ normal(0, 3)
Group-Level Effects:
~cond (Number of levels: 2)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 2.34 1.84 0.07 6.81 1.00 2565 1381
~subj (Number of levels: 8)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 29.28 1.50 26.35 32.31 1.00 3428 3109
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 779.28 15.21 749.95 808.98 1.00 1450 2043
patient1 -28.15 21.39 -70.74 14.90 1.00 1379 1888
condN -42.58 6.14 -54.73 -30.49 1.00 2989 2597
armND 1.93 3.11 -4.00 8.24 1.00 6042 2634
patient1:condN 2.31 6.34 -10.39 14.66 1.00 4705 3035
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 52.39 0.95 50.62 54.30 1.00 5083 2995
alpha 1.80 0.23 1.37 2.27 1.00 5431 3150
Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).In the Population-Level Effects section of Tab. 6 are now two estimates for the skewness parameter alpha (one for the controlled [alpha_Intercept] and one for the normal condition [alpha_condN]), while there is no entry for alpha in the Family Specific section anymore. All this due to our explicit modeling of alpha conditional on instruction condition.
pp_check(m5, nsamples = 100)
Fig. 11. Posterior predictive check. The thick blue line shows the distribution of the empirical data. The thin blue lines are one-hundred realizations of data generated from parameters estimated by the model.
conditional_effects(m5)
2.1.6 Group x Condition + ArmVI Model
(m6_form <-bf(transSC ~ patient * cond + arm +
(1 | subj) +
(1 | cond) +
(1 | arm)))transSC ~ patient * cond + arm + (1 | subj) + (1 | cond) + (1 | arm) (m6_prior <- m5_prior) prior class coef group resp dpar nlpar bound source
normal(0, 3) sd user
normal(0, 3) sigma user
normal(0, 2) alpha userif (MODEL) {
m6 <- brm(m6_form,
prior = m6_prior,
family = skew_normal(),
inits = "0",
data = drum_beats)
m6 <- add_criterion(m6,
"loo",
reloo = TRUE)
save(m6,
file = "m6.rda")
} else {
load("m6.rda")
}Tab. 7. Model summary.
(m6_summary <- summary(m6,
priors = TRUE)) Family: skew_normal
Links: mu = identity; sigma = identity; alpha = identity
Formula: transSC ~ patient * cond + arm + (1 | subj) + (1 | cond) + (1 | arm)
Data: drum_beats (Number of observations: 1102)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
Priors:
alpha ~ normal(0, 2)
Intercept ~ student_t(3, 735, 136.5)
sd ~ normal(0, 3)
sigma ~ normal(0, 3)
Group-Level Effects:
~arm (Number of levels: 2)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 2.42 1.81 0.12 6.72 1.00 2557 2041
~cond (Number of levels: 2)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 2.39 1.80 0.11 6.57 1.00 2758 1662
~subj (Number of levels: 8)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 29.30 1.54 26.43 32.44 1.00 3015 2588
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 779.65 15.49 750.27 810.64 1.00 1394 2089
patient1 -28.54 21.28 -70.50 12.05 1.00 1234 1812
condN -42.50 6.11 -54.84 -30.45 1.00 3088 2499
armND 2.02 5.30 -8.05 13.29 1.00 2642 1946
patient1:condN 2.41 6.36 -9.94 14.73 1.00 4722 3409
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 52.39 0.95 50.58 54.28 1.00 5124 3052
alpha 1.81 0.23 1.37 2.29 1.00 5198 3336
Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).In the Population-Level Effects section of Tab. 7 are now two estimates for the skewness parameter alpha (one for the controlled [alpha_Intercept] and one for the normal condition [alpha_condN]), while there is no entry for alpha in the Family Specific section anymore. All this due to our explicit modeling of alpha conditional on instruction condition.
pp_check(m6, nsamples = 100)
Fig. 12. Posterior predictive check. The thick blue line shows the distribution of the empirical data. The thin blue lines are one-hundred realizations of data generated from parameters estimated by the model.
wrap_plots(plot(conditional_effects(m6),
plot = FALSE)) +
plot_annotation(tag_levels = "A")
2.1.7 Group x Condition x ArmVI Model
(m6a_form <-bf(transSC ~ patient * cond * arm +
(1 | subj) +
(1 | cond) +
(1 | arm)))transSC ~ patient * cond * arm + (1 | subj) + (1 | cond) + (1 | arm) (m6a_prior <- m6_prior) prior class coef group resp dpar nlpar bound source
normal(0, 3) sd user
normal(0, 3) sigma user
normal(0, 2) alpha userif (MODEL) {
m6a <- brm(m6a_form,
prior = m6a_prior,
family = skew_normal(),
inits = "0",
data = drum_beats)
m6a <- add_criterion(m6a,
"loo",
reloo = TRUE)
save(m6a,
file = "m6a.rda")
} else {
load("m6a.rda")
}Tab. 8. Model summary.
(m6a_summary <- summary(m6a,
priors = TRUE)) Family: skew_normal
Links: mu = identity; sigma = identity; alpha = identity
Formula: transSC ~ patient * cond * arm + (1 | subj) + (1 | cond) + (1 | arm)
Data: drum_beats (Number of observations: 1102)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
Priors:
alpha ~ normal(0, 2)
Intercept ~ student_t(3, 735, 136.5)
sd ~ normal(0, 3)
sigma ~ normal(0, 3)
Group-Level Effects:
~arm (Number of levels: 2)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 2.41 1.82 0.08 6.76 1.00 2627 1431
~cond (Number of levels: 2)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 2.39 1.81 0.08 6.64 1.00 2860 1720
~subj (Number of levels: 8)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 29.40 1.53 26.51 32.57 1.00 3029 2344
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept 784.67 16.15 752.08 817.79 1.00 1483
patient1 -42.01 22.13 -85.49 0.98 1.00 1494
condN -42.41 7.62 -57.35 -27.23 1.00 2475
armND -7.47 7.37 -22.07 7.64 1.00 2409
patient1:condN 9.06 8.88 -8.67 26.35 1.00 2442
patient1:armND 26.58 8.72 9.34 43.94 1.00 2371
condN:armND -0.09 8.57 -16.94 16.33 1.00 2324
patient1:condN:armND -13.83 12.26 -37.64 10.59 1.00 2053
Tail_ESS
Intercept 1835
patient1 2077
condN 2221
armND 2811
patient1:condN 2681
patient1:armND 3043
condN:armND 2788
patient1:condN:armND 2859
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 52.18 0.93 50.37 54.05 1.00 4209 2815
alpha 1.69 0.24 1.21 2.17 1.00 4256 2965
Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).In the Population-Level Effects section of Tab. 8 are now two estimates for the skewness parameter alpha (one for the controlled [alpha_Intercept] and one for the normal condition [alpha_condN]), while there is no entry for alpha in the Family Specific section anymore. All this due to our explicit modeling of alpha conditional on instruction condition.
pp_check(m6a, nsamples = 100)
Fig. 13. Posterior predictive check. The thick blue line shows the distribution of the empirical data. The thin blue lines are one-hundred realizations of data generated from parameters estimated by the model.
wrap_plots(plot(conditional_effects(m6a),
plot = FALSE)) +
plot_annotation(tag_levels = "A")
conditions <- make_conditions(drum_beats,
vars = "arm")
conditional_effects(m6a,
effects = "patient:cond",
conditions = conditions)
2.1.8 Modeling Condition-specific Variation
As is apparent from Fig. 1, the controlled condition yielded broader density distributions in most descriptors. Hence, we model condition-dependent variation in the next model:
(m7_form <-bf(transSC ~ patient * cond * arm +
(1 | subj) +
(1 | cond) +
(1 | arm),
sigma ~ cond))transSC ~ patient * cond * arm + (1 | subj) + (1 | cond) + (1 | arm)
sigma ~ condThis model has now two outcomes, not just one, as the models before. Spread is estimated as sigma, and it varies conditional on instruction condition.
Here’s the non-standard (additional) prior:
(m7_prior <- c(set_prior("normal(0, 3)", class = "sd"),
set_prior("normal(0, 2)", class = "alpha"),
set_prior("normal(0, 2)", dpar = "sigma"))) prior class coef group resp dpar nlpar bound source
normal(0, 3) sd user
normal(0, 2) alpha user
normal(0, 2) b sigma userif (MODEL) {
m7 <- brm(m7_form,
prior = m7_prior,
family = skew_normal(),
inits = "0",
data = drum_beats)
m7 <- add_criterion(m7,
"loo",
reloo = TRUE)
save(m7,
file = "m7.rda")
} else {
load("m7.rda")
}Tab. 9. Model summary.
(m7_summary <- summary(m7,
priors = TRUE)) Family: skew_normal
Links: mu = identity; sigma = log; alpha = identity
Formula: transSC ~ patient * cond * arm + (1 | subj) + (1 | cond) + (1 | arm)
sigma ~ cond
Data: drum_beats (Number of observations: 1102)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
Priors:
alpha ~ normal(0, 2)
b_sigma ~ normal(0, 2)
Intercept ~ student_t(3, 735, 136.5)
Intercept_sigma ~ student_t(3, 0, 2.5)
sd ~ normal(0, 3)
Group-Level Effects:
~arm (Number of levels: 2)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 2.35 1.83 0.08 6.78 1.00 2801 1881
~cond (Number of levels: 2)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 2.35 1.79 0.09 6.66 1.00 2861 1790
~subj (Number of levels: 8)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 27.47 1.48 24.67 30.52 1.00 3279 2547
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept 788.48 15.51 758.86 818.79 1.00 1557
sigma_Intercept 4.31 0.04 4.24 4.38 1.00 4721
patient1 -40.96 21.12 -82.48 -0.43 1.00 1664
condN -45.95 7.92 -61.63 -30.64 1.00 2316
armND -9.51 9.12 -27.24 7.97 1.00 2081
patient1:condN 7.54 10.07 -11.46 27.47 1.00 2254
patient1:armND 23.88 11.56 1.26 46.68 1.00 2098
condN:armND 0.63 9.52 -17.93 19.30 1.00 2054
patient1:condN:armND -9.53 13.69 -35.98 17.18 1.00 2036
sigma_condN -0.47 0.06 -0.58 -0.36 1.00 4742
Tail_ESS
Intercept 2162
sigma_Intercept 2741
patient1 1948
condN 2890
armND 2631
patient1:condN 2928
patient1:armND 2658
condN:armND 2827
patient1:condN:armND 2641
sigma_condN 2671
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
alpha 2.38 0.30 1.83 3.00 1.00 4596 2823
Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).In the Links section of the summary (Tab. 9) we note that sigma is no longer modeled on the identity scale but on the \(\log_{2}\) scale. In the Population-Level Effects section of Tab. 9 we now find not only the estimates for the Intercept and condN but also two estimates for sigma (one for the controlled [sigma_Intercept] and one for the normal condition [sigma_condN]). Because of this, there is no entry for sigma in the Family Specific section anymore. All this due to our explicit modeling of sigma conditional on instruction condition.
pp_check(m7, nsamples = 100)
Fig. 14. Posterior predictive check. The thick blue line shows the distribution of the empirical data. The thin blue lines are one-hundred realizations of data generated from parameters estimated by the model.
p <- plot(conditional_effects(m7,
dpar = "mu"),
plot = FALSE)
p1 <- wrap_plots(p[[1]]+
p[[2]] +
p[[3]])
p <- plot(conditional_effects(m7,
dpar = "sigma"),
plot = FALSE)
p2 <- wrap_plots(p[[1]]+
p[[2]] +
p[[3]])
wrap_plots(p1, p2,
ncol = 1) +
plot_annotation(tag_levels = "A")
Fig. 15. Conditional plot.
In Fig. 15 we see a difference between the mean (‘mu’) estimated value for transSC depending on instruction condition (A), but no difference between dominant and non-dominant arm (B). Consequently, there is no interaction in (C).
The spread (sigma) of the estimated distributions also differs between conditions (D), but again not between arms (E). The latter is not surprising because we modeled sigma to vary conditional on condition, not arm.
2.1.9 Modeling Arm-specific Skewness
At least in some descriptors in Fig. 1 the skewness of the distribution seems to change depending on the arm. The skew normal distribution is a generalization of the Gaussian distribution, allowing for the additional shape-parameter skewness (asymmetry) to vary. Hence, we model this side-dependent skewness in the next model:
(m8_form <-bf(transSC ~ patient * cond * arm +
(1 | subj) +
(1 | cond) +
(1 | arm),
sigma ~ cond,
alpha ~ arm))transSC ~ patient * cond * arm + (1 | subj) + (1 | cond) + (1 | arm)
sigma ~ cond
alpha ~ arm(m8_prior <- c(set_prior("normal(0, 3)", class = "sd"),
set_prior("normal(0, 2)", dpar = "sigma"),
set_prior("normal(0, 2)", dpar = "alpha"))) prior class coef group resp dpar nlpar bound source
normal(0, 3) sd user
normal(0, 2) b sigma user
normal(0, 2) b alpha userif (MODEL) {
m8 <- brm(m8_form,
prior = m8_prior,
family = skew_normal(),
inits = "0",
data = drum_beats)
m8 <- add_criterion(m8,
"loo",
reloo = TRUE)
save(m8,
file = "m8.rda")
} else {
load("m8.rda")
}Tab. 10. Model summary.
(m8_summary <- summary(m8,
priors = TRUE)) Family: skew_normal
Links: mu = identity; sigma = log; alpha = identity
Formula: transSC ~ patient * cond * arm + (1 | subj) + (1 | cond) + (1 | arm)
sigma ~ cond
alpha ~ arm
Data: drum_beats (Number of observations: 1102)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
Priors:
b_alpha ~ normal(0, 2)
b_sigma ~ normal(0, 2)
Intercept ~ student_t(3, 735, 136.5)
Intercept_alpha ~ normal(0, 4)
Intercept_sigma ~ student_t(3, 0, 2.5)
sd ~ normal(0, 3)
Group-Level Effects:
~arm (Number of levels: 2)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 2.41 1.75 0.13 6.52 1.00 3444 1763
~cond (Number of levels: 2)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 2.39 1.88 0.06 6.86 1.00 2526 1619
~subj (Number of levels: 8)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 27.32 1.52 24.42 30.34 1.00 2796 2711
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept 788.57 15.53 758.00 819.87 1.00 1549
sigma_Intercept 4.34 0.04 4.26 4.41 1.00 4556
alpha_Intercept 1.73 0.38 1.00 2.45 1.00 4654
patient1 -44.64 21.25 -87.20 -3.22 1.00 1652
condN -44.62 8.57 -60.88 -27.20 1.00 2054
armND -9.99 9.34 -27.63 9.14 1.00 1745
patient1:condN 9.87 10.74 -11.39 31.30 1.00 2205
patient1:armND 27.84 12.01 3.26 50.58 1.00 2084
condN:armND -1.18 9.85 -21.34 17.80 1.00 1940
patient1:condN:armND -10.48 14.06 -38.01 17.74 1.00 1821
sigma_condN -0.51 0.06 -0.62 -0.40 1.00 4786
alpha_armND 1.71 0.53 0.71 2.80 1.00 4125
Tail_ESS
Intercept 1860
sigma_Intercept 3245
alpha_Intercept 2814
patient1 2043
condN 2373
armND 2146
patient1:condN 2608
patient1:armND 2325
condN:armND 2146
patient1:condN:armND 1918
sigma_condN 3009
alpha_armND 2658
Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).In the Population-Level Effects section of Tab. 10 are now two estimates for the skewness parameter alpha (one for the dominant [alpha_Intercept] and one for the non-dominant arm [alpha_armND]), while there is no entry for alpha in the Family Specific section anymore. All this due to our explicit modeling of alpha conditional on arm.
pp_check(m8, nsamples = 100)
Fig. 16. Posterior predictive check. The thick blue line shows the distribution of the empirical data. The thin blue lines are one-hundred realizations of data generated from parameters estimated by the model.
p <- plot(conditional_effects(m8,
dpar = "mu"),
plot = FALSE)
p1 <- wrap_plots(p[[1]] + xlab("") +
p[[2]] + xlab("") + ylab("") +
p[[3]] + xlab("") + ylab("") + theme(legend.position = "none"))
p <- plot(conditional_effects(m8,
dpar = "alpha"),
plot = FALSE)
p2 <- wrap_plots(p[[1]] + xlab("") +
p[[2]] + xlab("") + ylab("") +
p[[3]] + xlab("") + ylab("") + theme(legend.position = "none"))
p <- plot(conditional_effects(m8,
dpar = "sigma"),
plot = FALSE)
p3 <- wrap_plots(p[[1]] +
p[[2]] + ylab("") +
p[[3]] + ylab("") + theme(legend.position = "bottom"))
wrap_plots(p1, p2, p3,
ncol = 1) +
plot_annotation(tag_levels = "A")
Fig. 17. Conditional plots.
In Fig. 17 there is a clear difference in estimated mean (‘mu’) of \(transSC\) conditional on instruction (A), but not on arm (B). Consequently, there is no interaction between instruction and arm (C). On the other hand, there is a difference in estimated skewness (‘alpha’) between sides (E), but not instructions (D). (F) follows from that. ‘sigma’ shows a clear difference between instructions (G), not so for side (H), which is mirrored in (I).
2.1.10 Modeling Influence of Instruction and Arm on Shape
To be able to see whether skewness and spread are depending on both experimental manipulation and arm, the model formulas for ‘sigma’ and ‘alpha’ are updated accordingly:
(m9_form <-bf(transSC ~ patient * cond * arm +
(1 | subj) +
(1 | cond) +
(1 | arm),
sigma ~ cond + arm,
alpha ~ arm + cond))transSC ~ patient * cond * arm + (1 | subj) + (1 | cond) + (1 | arm)
sigma ~ cond + arm
alpha ~ arm + cond(m9_prior <- m8_prior) prior class coef group resp dpar nlpar bound source
normal(0, 3) sd user
normal(0, 2) b sigma user
normal(0, 2) b alpha userif (MODEL) {
m9 <- brm(m9_form,
prior = m9_prior,
family = skew_normal(),
inits = "0",
data = drum_beats)
m9 <- add_criterion(m9,
"loo",
reloo = TRUE)
save(m9,
file = "m9.rda")
} else {
load("m9.rda")
}Tab. 11. Model summary.
(m9_summary <- summary(m9,
priors = TRUE)) Family: skew_normal
Links: mu = identity; sigma = log; alpha = identity
Formula: transSC ~ patient * cond * arm + (1 | subj) + (1 | cond) + (1 | arm)
sigma ~ cond + arm
alpha ~ arm + cond
Data: drum_beats (Number of observations: 1102)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
Priors:
b_alpha ~ normal(0, 2)
b_sigma ~ normal(0, 2)
Intercept ~ student_t(3, 735, 136.5)
Intercept_alpha ~ normal(0, 4)
Intercept_sigma ~ student_t(3, 0, 2.5)
sd ~ normal(0, 3)
Group-Level Effects:
~arm (Number of levels: 2)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 2.43 1.85 0.09 6.87 1.00 2767 1986
~cond (Number of levels: 2)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 2.33 1.85 0.06 6.79 1.00 2605 1483
~subj (Number of levels: 8)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 27.27 1.56 24.40 30.38 1.00 3265 2317
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept 790.12 15.52 759.57 820.68 1.00 1742
sigma_Intercept 4.34 0.05 4.25 4.43 1.00 3828
alpha_Intercept 4.56 0.69 3.33 5.98 1.00 2811
patient1 -37.77 21.57 -80.96 5.12 1.00 1805
condN -46.73 7.65 -62.07 -31.63 1.00 2389
armND -15.95 8.85 -33.61 1.23 1.00 2545
patient1:condN 2.17 9.01 -15.34 19.64 1.00 2620
patient1:armND 32.98 10.44 13.01 53.36 1.00 2344
condN:armND 4.64 9.06 -13.35 22.42 1.00 2408
patient1:condN:armND -16.06 12.80 -41.41 9.03 1.00 2207
sigma_condN -0.62 0.05 -0.72 -0.51 1.00 4465
sigma_armND 0.06 0.05 -0.05 0.17 1.00 4040
alpha_armND 2.61 0.82 1.14 4.16 1.00 1376
alpha_condN -4.89 0.82 -6.49 -3.36 1.00 1558
Tail_ESS
Intercept 1784
sigma_Intercept 2963
alpha_Intercept 3012
patient1 2174
condN 2609
armND 2775
patient1:condN 2816
patient1:armND 2806
condN:armND 2826
patient1:condN:armND 2517
sigma_condN 3055
sigma_armND 3212
alpha_armND 3160
alpha_condN 2948
Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).In the Population-Level Effects section of Tab. 11 are now three estimates for the skewness parameter alpha, as well as three estimates for spread (‘sigma’).
pp_check(m9, nsamples = 100)
Fig. 18. Posterior predictive check. The thick blue line shows the distribution of the empirical data. The thin blue lines are one-hundred realizations of data generated from parameters estimated by the model.
p <- plot(conditional_effects(m9,
dpar = "mu"),
plot = FALSE)
p1 <- wrap_plots(p[[1]] + xlab("") +
p[[2]] + xlab("") + ylab("") +
p[[3]] + xlab("") + ylab("") +
p[[4]] + xlab("") + ylab("") + theme(legend.position = "none") +
p[[5]] + xlab("") + ylab("") + theme(legend.position = "none") +
p[[6]] + xlab("") + ylab("") + theme(legend.position = "none"))
p <- plot(conditional_effects(m9,
dpar = "alpha"),
plot = FALSE)
p2 <- wrap_plots(p[[1]] + xlab("") +
p[[2]] + xlab("") + ylab("") +
p[[3]] + xlab("") + ylab("") +
p[[4]] + xlab("") + ylab("") + theme(legend.position = "none") +
p[[5]] + xlab("") + ylab("") + theme(legend.position = "none") +
p[[6]] + xlab("") + ylab("") + theme(legend.position = "none"))
p <- plot(conditional_effects(m9,
dpar = "sigma"),
plot = FALSE)
p3 <- wrap_plots(p[[1]] +
p[[2]] + ylab("") +
p[[3]] + ylab("") +
p[[4]] + ylab("") + theme(legend.position = "bottom") +
p[[5]] + ylab("") + theme(legend.position = "bottom") +
p[[6]] + ylab("") + theme(legend.position = "bottom"))
wrap_plots(p1, p2, p3) +
plot_annotation(tag_levels = "A")
Fig. 19. Conditional plots.
In Fig. 19 there is a clear difference in estimated mean (‘mu’) of \(transSC\) conditional on instruction (A), but not on arm (B). Consequently, there is no interaction between instruction and arm (C). There is a difference in estimated skewness (‘alpha’) both between instructions (D) and sides (E), but with opposing trends. Nevertheless, there is no interaction between them (F). ‘sigma’ shows a clear difference between instructions (G) as well as for side (H), but albeit their opposing trends, there is no interaction (I).
2.1.11 Model Comparison
We compare models by their estimated log-posterior density (elpd). The smaller this value the better a model predicts the data, despite the penalty for additional covariates. When the difference between two models is more than two SE apart they are considered to be ‘different enough’ to warrant acceptance of one over the other despite possibly smaller parsimony.
loo_compare(m0, m2, m3, m5, m6, m6a, m7, m8, m9) elpd_diff se_diff
m9 0.0 0.0
m8 -38.1 8.0
m7 -42.2 8.9
m6a -80.4 16.0
m2 -80.7 16.9
m3 -80.9 16.9
m6 -82.6 17.0
m5 -82.6 17.0
m0 -119.6 15.0 The model with the smallest elpd is m9, which corresponds to the last model.
2.1.12 Conclusion
Each row of the plot in Fig. 19 represents the results of one of the three outcomes, each column shows the estimates for one condition/manipulation. The dots in these plots are the means of the posterior distributions of the respective estimate, and can be interpreted similarly to the point estimates in traditional frequentist statistics. The error bars represent 95% credible intervals, which are interpreted as comprising the true value with 95% probability, given the model and the data.
From the first row of Fig. 19 we can infer that instructions had an effect on the mean, while arm had not. The second row suggests that instruction had a larger, and opposing, effect on skewness compared to arm arm, and finally, the third row claims the same for spread.
On a different level of understanding these results this means that transSC is lower in the normal strokes, with almost no skewness (alpha is close to zero) and lower spread than in controlled strokes.
Based on the above we can infer that instruction causes the most pronounced differences in transSC, while the playing arm only seems to invoke minor to no differences.
2.2 Adding the Response transFlat
In the following, the descriptor transFlat is added as an outcome to our model. Technically speaking, we would need to re-start our modeling attempt with the Null model and work our way up to the ‘best’ model as we have done with transSC. This is necessary because we cannot just assume that with an additional outcome (which will actually become three additional outcomes: mu, alpha and sigma) will be best explained by the same model as was with just one outcome variable.
But we take a short cut here and just add transFlat to the best model of transSC.
To be able to see whether skewness and spread are depending on both experimental manipulation and arm, the model formulas for ‘sigma’ and ‘alpha’ are updated accordingly:
(m10_form <-bf(mvbind(transSC, transFlat) ~ cond + arm +
(1|p|cond) +
(1|q|arm) +
(1|r|subj),
sigma ~ cond + arm,
alpha ~ arm + cond))transSC ~ cond + arm + (1 | p | cond) + (1 | q | arm) + (1 | r | subj)
sigma ~ cond + arm
alpha ~ arm + cond
transFlat ~ cond + arm + (1 | p | cond) + (1 | q | arm) + (1 | r | subj)
sigma ~ cond + arm
alpha ~ arm + condThe prior probability on the individual intercepts needs to reflect that there are now two response variables:
(m10_prior <- c(
set_prior("lkj(2)", class = "cor"),
set_prior("normal(0, 0.1)", class = "sd", resp = "transFlat"),
set_prior("normal(0, 2)", class = "b", resp = "transFlat", dpar = "alpha"),
set_prior("normal(0, 2)", class = "Intercept", resp = "transFlat", dpar = "alpha"),
set_prior("normal(0, 1)", class = "Intercept", resp = "transFlat", dpar = "sigma"),
set_prior("normal(0, 1)", class = "b", resp = "transFlat", dpar = "sigma"),
set_prior("normal(0, 0.1)", class = "sd", resp = "transSC"),
set_prior("normal(0, 2)", class = "b", resp = "transSC", dpar = "alpha"),
set_prior("normal(0, 2)", class = "Intercept", resp = "transSC", dpar = "alpha"),
set_prior("normal(0, 2)", class = "Intercept", resp = "transSC", dpar = "sigma"),
set_prior("normal(0, 2)", class = "b", resp = "transSC", dpar = "sigma")
)) prior class coef group resp dpar nlpar bound source
lkj(2) cor user
normal(0, 0.1) sd transFlat user
normal(0, 2) b transFlat alpha user
normal(0, 2) Intercept transFlat alpha user
normal(0, 1) Intercept transFlat sigma user
normal(0, 1) b transFlat sigma user
normal(0, 0.1) sd transSC user
normal(0, 2) b transSC alpha user
normal(0, 2) Intercept transSC alpha user
normal(0, 2) Intercept transSC sigma user
normal(0, 2) b transSC sigma userif (MODEL) {
t <- Sys.time()
m10 <- brm(m10_form,
prior = m10_prior,
family = skew_normal(),
data = drum_beats,
inits = "0")
print(Sys.time() - t)
m10 <- add_criterion(m10,
"loo",
reloo = TRUE)
save(m10,
file = "m10.rda")
} else {
# load("m10.rda")
}NULLTab. 12. Model summary.
# (m10_summary <- summary(m10,
# priors = TRUE))In the Population-Level Effects section of Tab. 12 are now three estimates for the skewness parameter alpha, as well as three estimates for spread (‘sigma’).
# pp_check(m10, resp = "transSC", nsamples = 100)
# pp_check(m10, resp = "transFlat", nsamples = 100)Fig. 20. Posterior predictive check. The thick blue line shows the distribution of the empirical data. The thin blue lines are one-hundred realizations of data generated from parameters estimated by the model.
# p <- plot(conditional_effects(m10,
# dpar = "mu"),
# plot = FALSE)
# p1 <- wrap_plots(p[[1]] + xlab("") +
# p[[2]] + xlab("") + ylab("") +
# p[[3]] + xlab("") + ylab("") + theme(legend.position = "none"))
# p <- plot(conditional_effects(m10,
# dpar = "alpha"),
# plot = FALSE)
# p2 <- wrap_plots(p[[1]] + xlab("") +
# p[[2]] + xlab("") + ylab("") +
# p[[3]] + xlab("") + ylab("") + theme(legend.position = "none"))
# p <- plot(conditional_effects(m10,
# dpar = "sigma"),
# plot = FALSE)
# p3 <- wrap_plots(p[[1]] +
# p[[2]] + ylab("") +
# p[[3]] + ylab("") + theme(legend.position = "bottom"))
# wrap_plots(p1, p2, p3,
# ncol = 1) +
# plot_annotation(tag_levels = "A")Fig. 21. Conditional plots.
In Fig. 21 there is a clear difference in estimated mean (‘mu’) of \(transSC\) conditional on instruction (A), but not on arm (B). Consequently, there is no interaction between instruction and arm (C). There is a difference in estimated skewness (‘alpha’) both between instructions (D) and sides (E), but with opposing trends. Nevertheless, there is no interaction between them (F). ‘sigma’ shows a clear difference between instructions (G) as well as for side (H), but albeit their opposing trends, there is no interaction (I).
2.3 transCrest
3 System Setup
Data wrangling and analyses were carried out with the statistical package R (R version 4.0.2 (2020-06-22); R Core Team 2020). Bayesian modeling was done with the package brms (Bürkner 2017, 2018) which uses the probabilistic language Stan as back end (Carpenter et al. 2017). Plots were done with the packages bayesplot (Gabry and Mahr 2020) and ggplot2 (???).
References
Bates, D. M. 2010. Lme4: Mixed-Effects Modeling with R. Springer. http://lme4.r-forge.r-project.org/lMMwR/lrgprt.pdf.
Bürkner, Paul-Christian. 2017. “Brms: An R Package for Bayesian Multilevel Models Using Stan.” Journal of Statistical Software 80 (1). https://doi.org/10.18637/jss.v080.i01.
———. 2018. “Advanced Bayesian Multilevel Modeling with the R Package Brms.” The R Journal 10 (1).
Carpenter, Bob, Andrew Gelman, Matthew Hoffman, Daniel Lee, Ben Goodrich, Michael Betancourt, Marcus Brubaker, Jiqiang Guo, Peter Li, and Allen Riddell. 2017. “Stan: A Probabilistic Programming Language.” Journal of Statistical Software, Articles 76 (1): 1–32. https://doi.org/10.18637/jss.v076.i01.
Câmara, Guilherme Schmidt, Kristian Nymoen, Olivier Lartillot, and Anne Danielsen. 2020. “Effects of Instructed Timing on Electric Guitar and Bass Sound in Groove Performance.” The Journal of the Acoustical Society of America 147 (2): 1028–41. https://doi.org/10.1121/10.0000724.
Danielsen, Anne, Carl Haakon Waadeland, Henrik G. Sundt, and Maria A. G. Witek. 2015. “Effects of Instructed Timing and Tempo on Snare Drum Sound in Drum Kit Performance.” The Journal of the Acoustical Society of America 138 (4): 2301–16. https://doi.org/10.1121/1.4930950.
Gabry, Jonah, and Tristan Mahr. 2020. “Bayesplot: Plotting for Bayesian Models.” https://mc-stan.org/bayesplot.
R Core Team. 2020. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.
Tibshirani, Robert. 1996. “Regression Shrinkage and Selection via the Lasso.” Journal of the Royal Statistical Society: Series B (Methodological) 58 (1): 267–88. https://doi.org/10.1111/j.2517-6161.1996.tb02080.x.
Wilkinson, G. N., and C. E. Rogers. 1973. “Symbolic Description of Factorial Models for Analysis of Variance.” Applied Statistics 22 (3): 392–99. https://doi.org/10.2307/2346786.